%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Langevin equation without approximation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Langevin

Constant;
kap=1;
N  =1e6;									Ns =200;
U0 =-.1e-6*kap;							T  =650e-6;
eta=100000*kap;							eps=.33*kap;
gFr=10*kap;								Dif=(4*pi/797e-9/2/pi/22.5e3)^2*kB*T/85/uA/gFr;
dt=2e-3/kap;

t =   0;
ap= eta;									r=rand(1,Ns)*15;
am=   0;									p=sqrt(50*T/eps)*randn(1,Ns);

figure(3);
fig=figure;
set(3,'DoubleBuffer','on');%'Position',[10 10 1010 680],
set(gca,'NextPlot','replace','Visible','off');
%mov = avifile('C:\Dokumente und Einstellungen\Philippe\Eigene Dateien\Vortraege\Carl\TalkLangevin_1.avi')

skip=50;									leng=10000;
for jj=1:leng
	if((jj/skip)==floor(jj/skip))
		Am(jj/skip)=am*am';				rCoM(jj/skip)=1/Ns*sum(r);
		Ab(jj/skip)=2*ap*am'+2*ap'*am;	bnch(jj/skip)=abs(1/Ns*sum(exp(2*i*r)));
		tt(jj/skip)=t;					Fm(jj/skip)=am;
jj/skip
		LangevinFigure(ap,am,r,Am,Ab,rCoM,bnch,tt);
%		eval(['print -dbitmap D:\Fe',num2str(jj/skip+100)])
%        F  =getframe(gcf);
%        mov=addframe(mov,F);
	end
	te = sqrt(2*Dif*dt)*randn(1,Ns);

	t  = t  + dt;
	b0 = -kap*am - i*U0*N/Ns*sum(exp(2*i*r ))*ap + eta*(t<0);
	p0 = p;
	f0 = 4*i*eps*U0*ap*(am*exp(-2*i*r )-am'*exp(2*i*r )) - gFr*p*(t>5) ;

	a1 = am + dt*b0;
	r1 = r  + dt*p0*eps;
	p1 = p  + dt*f0 + te;

	b1 = -kap*a1 - i*U0*N/Ns*sum(exp(2*i*r1))*ap + eta*(t<0);
	p1 = p1;
	f1 = 4*i*eps*U0*ap*(a1*exp(-2*i*r1)-a1'*exp(2*i*r1)) - gFr*p1*(t>5);

	am = am + dt*.5*(b0+b1);
	r  = r  + dt*.5*(p0+p1);
	p  = p  + dt*.5*(f0+f1) + te;
end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subroutine.

function LangevinFigure(ap,am,r,Am,Ab,rCoM,bnch,tt)

kap=1;%2*pi*22.5e3;
x =[0:.005:5]*pi;
figure(3);		clf;
subplot(221),plot(tt/kap,Am,'LineWidth',1.5)
set(gca,'xLim',[0 20],'yLim',[0 1e7],'XTickLabel',' ','YTickLabel',' ','LineWidth',1.3);
xlabel('\fontsize{18}{\ittime}')
ylabel('\fontsize{18}{\itprobe power}')
subplot(222),plot(tt/kap,bnch,'LineWidth',1.5)
set(gca,'xLim',[0 20],'yLim',[0 1],'XTickLabel',' ','FontSize',14,'LineWidth',.31);
xlabel('\fontsize{18}{\ittime}')
ylabel('\fontsize{18}{\itbunching}')
subplot(224),plot(tt/kap,Ab,'k','LineWidth',1.5)
set(gca,'xLim',[0 20],'yLim',[-2e9 2e9],'XTickLabel',' ','YTickLabel',' ','LineWidth',1.3);
xlabel('\fontsize{18}{\ittime}')
ylabel('\fontsize{18}{\itbeat note}')
EE=abs(ap*exp(i*x)+am*exp(-i*x)).^2;
EP=abs(ap*exp(i*r)+am*exp(-i*r)).^2;
subplot(223),plot(mod(r/pi,5),-EP+6e9,'.b',x/pi,-EE+6e9,'r','MarkerSize',15,'LineWidth',1.5)
set(gca,'xLim',[0 5],'yLim',[-6e9 0],'XTickLabel',' ','YTickLabel',' ','LineWidth',1.3);
xlabel('\fontsize{18}{\itatom distribution}')
text(1.5,-1.5e9,'\fontsize{18}pump \rightarrow')
pause(.01)

end
